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Over the past eight hundred thousand years, glacial—interglacial cycles oscillated witha 
period of one hundred thousand years (‘100k world”). Ice core and ocean sediment data 
have shown that atmospheric carbon dioxide, Antarctic temperature, deep ocean 
temperature, and global ice volume correlated strongly with each other in the 100k 
world” ®. Between about 2.8 and 1.2 million years ago, glacial cycles were smaller in 
magnitude and shorter in duration (40k world”). Proxy data from deep-sea sediments 
suggest that the variability of atmospheric carbon dioxide in the 40k world was also 
lower than in the 100k world®°, but we do not have direct observations of atmospheric 
greenhouse gases from this period. Here we report the recovery of stratigraphically 
discontinuous ice more than two million years old from the Allan Hills Blue Ice Area, East 
Antarctica. Concentrations of carbon dioxide and methane inice core samples older 
than two million years have been altered by respiration, but some younger samples are 
pristine. The recovered ice cores extend direct observations of atmospheric carbon 
dioxide, methane, and Antarctic temperature (based on the deuterium/hydrogen 
isotope ratio 5D,,.., a proxy for regional temperature) into the 40k world. All climate 
properties before eight hundred thousand years ago fall within the envelope of 
observations from continuous deep Antarctic ice cores that characterize the 100k 
world. However, the lowest measured carbon dioxide and methane concentrations and 
Antarctic temperature in the 40k world are well above glacial values from the past eight 
hundred thousand years. Our results confirm that the amplitudes of glacial—interglacial 
variations in atmospheric greenhouse gases and Antarctic climate were reduced in the 
40k world, and that the transition from the 40k to the 100k world was accompanied by a 


decline in minimum carbon dioxide concentrations during glacial maxima. 


Earth has been cooling, and ice sheets expanding, over approximately 
the past 52 million years (Myr)". Superimposed on this cooling are peri- 
odic changes in the Earth’s climate system driven by variations in the 
eccentricity (with periods of 400 and 100 kyr) and precession (23 and 
19 kyr) of the Earth’s orbit around the Sun, and the tilt of the spin axis 
(about 41 kyr). From around 2.8 to 1.2 Myr ago (Ma), the Earth’s climate 
system oscillated between glacial and interglacial states with a period 
of about 40 kyr (the 40k world’). Between 1.2 and 0.8 Ma, an interval 
known as the ‘mid-Pleistocene transition’ (MPT), the period of glacial 
cycles lengthened to about 100 kyr and glacial periods became colder”, 
for reasons that are poorly understood. The post-MPT glacial cycles are 
characterized by a quasi-100-kyr period (the 100k world’). 

Studies of stratigraphically continuous ice cores have shown that 
atmospheric CO, is directly linked to Antarctic and global temperature 
over the last 800 kyr*ć. The coupling of the Earth’s climate and carbon 
cyclein earlier times, however, has not been conclusively demonstrated. 
Recent boron-based reconstructions have suggested that glacial CO, 


declined across the MPT, but interglacial CO, remained stable*°. These 
studies suggest that the range of CO, was reduced in the 40k world, as 
was the amplitude of glacial cycles. On the other hand, CO, estimates 
based on foraminiferal 85°C predict that CO, varied between 300 and 180 
parts per million (ppm) in the 40k world®”, hinting at a fundamentally 
different relationship between CO, and temperature. As both of these 
techniques rely on indirect measurements of atmospheric CO,, they 
have lower accuracy and precision (typically more than + 20 ppm”; 
2s.d.(o)) than can beachieved by direct measurements inice core samples 
(+2 ppm; 20). Accurate and precise CO, measurements from ice cores 
that pre-date the MPT are therefore needed. 


Old ice records from the Allan Hills 


Million-year-old ice was recently discovered” in shallow ice cores drilled 
in the Allan Hills Blue Ice Area (BIA), Antarctica (—76.73° N, 159.36° E; 
Extended Data Figs. 1, 2). In the Allan Hills, bedrock topography and 
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Fig. 1| Age-depth profile of Allan Hills ice cores. a, ALHIC1502 data in blue; 

b, ALHIC1503 data in red. c, *°Ar,,, ages plotted against distance above bedrock 
with the same colour coding asin a, b. Data in circles were measured at Princeton 
University”. New measurements made at Scripps Institution of Oceanography are 
plotted as squares. Error bars represent the external reproducibility (10;110 kyr 
and 220 kyr for samples measured at Scripps and Princeton, respectively) ofthe 
measurement or 10% of the sample age, whichever is greater. Shading highlights 
sections in which ice that is more than 1 Myr old is present. Sections affected by 
respiration are marked by the vertical arrows (blue, ALHIC1502; red, ALHIC1503). 


strong katabatic winds lead to the exhumation of old ice from depth to 
the surface”. In 2015-2016, two additional cores were drilled. The first 
extended the drilling to a depth of 147 m at the BIT-58 borehole (named 
ALHIC1503 hereafter), where approximately 1-Myr-old ice had previously 
been discovered. The second (ALHIC1502), drilled 56 mto the southwest, 
reached bedrock at 191 m. These ice cores are dated by measuring the Ar 
deficit in trapped air relative to the modern atmosphere” (Ar m), and 
represent an archive of the Earth’s climate extending well into the 40k 
world. A nearby blue-ice core (Allan Hills Site 27; S27) gives a continuous 
record from about 100 to 250 ka, and serves as a reference section for 
comparison with ALHIC1502 and ALHIC1503". 

Age-depth profiles of the ALHIC1502 and ALHIC1503 cores (Fig. 1) 
show that both sites contain 300-500-kyr-old ice, overlying a basal unit 
confined to approximately the bottom 30 m, with ages that range from 
more than 1 to 2.7 Myr. At ALHIC1503, the deepest sample lies within 
about 1 m of bedrock and has the oldest *°Ar ım age, 2.7 + 0.3 (10) Myr. 
However, stratigraphic disturbance in both ALHIC1502 and ALHIC1503 
is evident from abrupt age discontinuities (Fig. 1), accompanied by large 
swings in 80 of O, (6°0,,,) and 5D,,. values”. In light of these compli- 
cations, we interpret the paleoclimate records from the Allan Hills BIA 
as discrete ‘snapshots’—discontinuous but accurate portraits—of the 
climate state, rather than continuous time series. 


Climate properties in the old ice 


Paleoclimate reconstructions from stratigraphically discontinuous ice 
sections face several challenges. First, there is the potential for sampling 
bias inthe preserved climate states due to possible differences in accu- 
mulation rate and/or ice rheology. Second, diffusion in million-year-old 
ice may lead to smoothing of paleoclimate properties or proxies’””°, 
though this effect is likely to be minor in Allan Hills ice (see Methods). 
Third, discrete sampling might not capture the complete glacial- 
interglacial range of properties of interest. Overall, we estimate that 
76 +14% (20) of the total atmospheric CO, variability in the 40k world is 
recovered by our CO, samples. This estimate is based upon our deduc- 
tions that 90% of the true glacial-interglacial climate variability is 
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Fig. 2| Climate properties over the past 2.9 Myr documented in ice core and 
benthic foram records. a, The continuous 800-kyr 5D,.. record from the 
European Project for Ice Coring in Antarctica, Dome C (EDC) ice core* (black 
line); the continuous S27 6D;,. record covering 120-250 ka" (red line); and the 
discrete Allan Hills 5D,.. ice samples (red circles). b, The 800-kyr ice core CH, 
record” (orange line) and the binned Allan Hills CH, data (orange circle). c, The 
800-kyr ice core CO, record*>*”* (purple line) and the binned Allan Hills CO, data 
(purple circles). Note that there are no reliable CO, and CH, analyses inthe 
2.7-Ma bin (see Methods). d, The 800-kyr ice core 8°O,m record”? (brown line) 
and binned values of Allan Hills 6*80,,,, samples (brown circles). All 880, Values 
are normalized to the modern atmosphere. e, The globally distributed benthic 
oxygen isotope stack over the Plio-Pleistocene (LRO4)” showing decreased 
glacial-interglacial variability and less extreme glacials before 800 ka. 


preserved in the Allan Hills ice, and that 84 + 15% (20) of CO, variability 
within the ice is captured by discrete samples, assuming a linear relation- 
ship between CO, and benthic 50 (see Methods). 

Measured climate proxies and ice core properties across time (Fig. 2) 
include proxy records of regional temperature (6D,,..), CH,, CO,, and global 
O, fractionation (8°0O,,,,)”. Samples older than 800 kyr are assigned the 
age of the nearest *°Ar,.m-dated sample and then binned into two age 
groups: the MPT (800-1,200 kyr), and the 40k world (1.2-2.7 Myr). Within 
the 40k world bin, samples are further subdivided into three groups with 
average ages of 1.5, 2.0, and 2.7 Myr (see Methods and Fig. 2). 

Improbably high CO, concentrations or 6°C values of CO, indicat- 
ing contributions from respiratory CO, were present in all 2.7-Myr-old 
samples, common among 2.0-Myr-old samples, but absent from 1.5-Myr- 
old samples (see Methods). We chose -7%o as the cut-off value for °C 
and rejected all CO, samples below the shallowest °C sample with an 
isotopic value of lighter than —7%o. The application of these criteria 
excluded 21 CO, samples from our analysis. CH, and 8O,.m samples 
were also rejected when 6°C or CO, measurements indicated that 
these properties were compromised. We continued to plot 6D,,, val- 
ues of the ice affected by respiration in Fig. 2 in the 2.0 and 2.7 Myr age 
groups, but excluded 8D,,, data from anomalous CO, samples in the 
40k world. 
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Fig. 3 | Comparison of blue-ice CO, record and boron-based CO, 
reconstructions during and before the MPT. a, LRO4 benthic oxygen isotope 
stack’. b, Distribution of MPT CO, (left, orange) and pre-MPT CO, (right, red) 
data observed in Allan Hills ice. c, Comparison between ice core CO, (orange and 
red circles) and 6"B-based CO, reconstructions (purple and blue solid lines; 


Ice dating back to about 100-250 ka, sampled at a nearby site’®, records 
a6D,,, range of -280 to -360%o in the 100k world. 5D,,, values ranged 
from —284 to -316% in 40k-world ice and from -288 to -332% in MPT ice. 
Maximum values of 6D,,, were similar during all three time periods, lying 
within an 8% range. However, minimum values were higher in the MPT 
and the 40k world by at least 28%o. This implies that interglacial periods 
during the 40k world in Antarctica were not significantly warmer than 
those inthe 100k world, but that glacial maxima were lower after the MPT. 

CO, concentrations ranged from 221 to 277 ppm (n= 23)" in MPT ice, 
and from 214 to 279 ppm (n = 37) in ice from the 40k world. Given that 
76 +14% (20) of total CO, variability was captured by 37 discrete samples, 
and assuming that glacial and interglacial extremes are equally absent, 
the expected true range of 40k world CO, is 8613 ppm (20; Fig. 3). Our 
best estimate of 40k-world CO, concentrations is therefore 204-289 ppm. 
These ranges indicate that interglacial CO, concentrations in the 40k and 
100k worlds were similar, whereas glacial CO, concentrations are likely 
to have been 24 ppm higher in the 40k world than in the 100k world. By 
comparison, reconstructions of atmospheric CO, from individual boron 
isotope measurements in foraminifera ranged from 190 to 320 ppm 
(Fig. 3), with the average extreme values being 234 and 277 ppm”. 


Implications for climate evolution 


Our findings appear to be inconsistent with hypotheses that attrib- 
ute the transition into the 100k world to a long-term decline in both 
interglacial and glacial atmospheric CO,””. Our data instead support 
hypotheses that link the MPT to greater ice sheet size and CO, drawdown 
during glacial maxima. Those hypotheses include, but are not limited to, 
enhanced dust delivery to the southern ocean” and changes in global 
ocean circulation” that resulted in an additional drawdown of atmos- 
pheric CO,, thereby enhancing the growth of continental ice sheets and 
leading to the emergence of 100-kyr glacial cycles. 

Compared to CO, concentrations and 6D,,., CH, concentrations and 
&'8O.,.m Values in ice from the MPT and 40k world exhibited even less 


dashed lines represent the 95% confidence intervals)”. Shading represents the 
expanded range estimate, given the possibility that discrete Allan Hills CO, 
samples may not fully capture the true range of CO, variability. Black dashed 
lines represent the glacial-interglacial range of CO, inthe 100k world’. 


glacial-interglacial variability. CH, concentrations ranged from 405 
to 569 parts per billion (ppb) in MPT ice (n = 31) and 457 to 592 ppb in 
40k-world ice (n = 16). These ranges are only about 40% of the range 
in the 100k world (approximately 350 to 800 ppb)”. Undersampling 
may play a role in this reduction, as CH, maxima in the 100k world have 
short durations”®. In addition, discrete CH, measurements require 
more ice than CO, measurements (60 g versus 10 g). They therefore 
average over a longer interval of time than CO, samples. &'°0,,,, val- 
ues ranged from 0.01% to 0.92% in MPT ice (n = 63) and from 0.00%o 
to 0.33% in 40k-world ice (n = 29). In comparison, the range is from 
—0.30%o to +1.48%o in ice from the 100k world”. The 6'%0,,,, values in 
ice from the 40k world are remarkable in that they lack values higher 
than +0.5%o and the total range of &'%0,,,, Values is only about 20% of 
that observed in ice from the 100k world. Smaller changes in seawater 
880 values (Fig. 2) help to explain the absence of 67%0,,,, values higher 
than +0.5%o. However, additional mechanisms are likely to be needed 
to explain the lack of precession-modulated variability in 5°O,, associ- 
ated with biosphere productivity and changes in the low-latitude water 
cycle”. 

Plots of atmospheric CO, against 5D,,. and CH, against 5D,,, (Fig. 4) 
show that samples fromthe MPT and 40k world fall within the envelope 
of the 100k world but exhibit only a fraction of its range. Samples from 
the MPT and 40k world do not have the isotopically light 6D,,. and low 
CO, values characteristic of glacial samples in the 100k world. Both CO, 
and CH, are positively correlated with 5D,,., in ice from the 100k world 
(S27). By contrast, in the MPT and 40k world, only 5D;,. and CO, are 
significantly correlated (P= 0.002 and P= 0.003, respectively). For all 
three sample sets, the slopes of the correlation between CO, and 5D,,.. 
are statistically indistinguishable (Fig. 4), in agreement with previous 
reconstructions from boron-based proxies and benthic 80 values. 
Together, these observations suggest a consistent and persistent cou- 
pling between Antarctic climate and atmospheric CO, throughout the 
Pleistocene®”’, and that glacial cycles in the 40k world are truncated 
versions of glacial cycles in the 100k world. 
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Fig. 4 | Cross-correlations between climate properties. a, CO, versus 8D;ce. 

b, CH, versus 6D,,.. Data points are colour-coded by age: Black, Vostok CO, or 
CH, data? projected onto coeval Allan Hills 5D,,, from S27; orange, Allan Hills 
MPT data; red, Allan Hills 40k world data. Shading represents the range of 6Dice 
we observed in the 40k world. 5D,.. measurements in samples measured for CO, 
span 94% of the total range of 5D,.. observed in 40k world samples, and therefore 
the CO, samples represent nearly the full range of climate variability recorded in 
the cores. There are significant correlations between CO, and 6D,,. in all three 
age units, and their regression slopes (ppm by %v; solid lines) are statistically 
indistinguishable: 1.33 + 0.17 (20; r= 0.78), 1.14 + 0.68 (r=0.61), and 0.90 + 0.56 
(r=0.48) for 120-250 ka, MPT, and pre-MPT age units, respectively. For 
comparison, a significant (P< 0.05) correlation between CH, and 6D,,, exists 
only at S27 (3.56 + 0.65; r= 0.61). 


Conclusions 


Shallow ice cores drilled in the Allan Hills BIA provide a direct observation 
of the variability in atmospheric CO, during the 40k world, confirming 
previous findings that minimum CO, concentrations declined after the 
MPT®”. These results complement plans to drill a continuous ice core 
record back to around 1.5 Ma”. These snapshots of the Earth’s climate 
system from shallow ice cores in BIAs provide an incomplete picture, 
especially with regards to the dynamics. Nevertheless, the oldest ice will 
probably be found in discontinuous sections. Our work demonstrates 
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that BIAs can be exploited to extend climate records, including atmos- 
pheric CO, concentrations, well back into the early Pleistocene and 
perhaps even the Pliocene. 
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Methods 


Sample location and description 

Allan Hills Blue Ice Area (BIA). The Allan Hills BIA is located -100 km 
to the northwest of the McMurdo Dry Valleys. Drilling sites are located 
to the west of the Allan Hills nunatak in the Main Ice Field (MIF) of the 
Allan Hills BIAs (Extended Data Fig. 1). Here, net ablation leads to the 
exposure of ancient glacial ice at the surface of an ice sheet”. Elevation, 
bedrock topography, and ice velocities in the MIF have been previously 
documented’. Site meteorology and surface snow properties are 
discussed in more detail by Dadic et al.* and ice and gas properties 
described in Higgins et al.”. 


Site ALHIC1503 (previously named BIT-58). Site ALHIC1503 
(76.73243° S, 159.3562° E) is located near a crest of anorthwest-southeast 
trendingice ridge off the MIF. The precise direction and the velocity of ice 
flow are not known at this location. Measurements from the nearby ice 
field are consistent with primary flow to the east or northeast. Complex 
dust bands are visible from high-resolution satellite imagery (Extended 
Data Fig. 1, left), implying a complex history of glacial flow. Typical ice 
flow velocities near ALHIC1503 are low (0.015 m yr)”. 

Bedrock topography, determined using ground penetrating radar 
(GPR), indicates that ALHIC1503 is located ona steep (~45°) slope lead- 
ing to a local bedrock high with an overlying ice thickness of -150 m 
(Extended Data Fig. 2). In the 2015-2016 field season, 23 m of ice was 
drilled, extending the earlier 124-m-long ice core” to the bedrock at 
147 m below the ice surface. 


Site ALHIC1502. Site ALHIC1502(76.73286°S, 159.35507° E) islocated 56m 
upflow from ALHIC1503. The ice thickness is -204 m, consistent witha 
steep bedrock slope (Extended Data Fig. 2). Drilling at site ALHIC1502 
recovered 197 m ofice. 


Site 27 (S27). S27 (76.70° S, 159.31° E) is located along the main ice flow 
line of the Allan Hills region. Two hundred and twenty-four metres of ice 
was retrieved, representing ~70% of the estimated column thickness. 
Spaulding et al. provide a comprehensive suite of stable water isotope 
and gas analyses for S27. The S27 ice chronology has been established by 
matching 6D,,, features to those of the Vostok ice core. This timescale 
was validated by correlating variations in 50,,,, (6'°O of paleoatmos- 
pheric O,) into the record of 8°0,,,, in Vostok ice cores”. 


Analytical procedures 

4° Ar stm and SXe/Kr. *°Ar is produced in the solid earth by the radioactive 
decay of *°K. As *°Ar slowly leaks into the atmosphere, its concentration 
increases with time. Ar and *Ar, on the other hand, are stable, primordial, 
and non-radiogenic, so their atmospheric concentrations are treated as 
constant. The ratio of “Ar/*Ar thus rises towards the future and decreases 
towards the past. The term of merit hereis the gravitationally corrected paleo- 
atmospheric *°Ar/*Ar ratio: Af gm = 5° Ar/*Ar—88Ar/*Ar. 

The latter term corrects for the gravitational fractionation in the firn 
column®. Studies of *°Ar m as a function of age in the Dome C and Vostok 
ice cores constrained the rate of change to be 0.066 + 0.006% Myr”, 
allowing us to date the air trapped in the ice without the prerequisite 
for stratigraphic continuity”. 

A potential complication of *°Ar,,,, dating comes from the in situ 
production of *°Ar from the bedrock and/or from dust, a tiny 
amount of which is present in polar ice cores. The deepest sample 
at ALHIC1502 has a much younger age (0.8 Myr) than the overlying 
ice (2.2 Myr). This inverted age—depth relationship may result either 
from folding or from radiogenic input of “Ar from the“’K in the underlying 
bedrock, a phenomenon observed for basal ice at GISP2 in Greenland”. 
The “°K from dust in the ice would produce a negligible amount of *°Ar. 

Ar isotope analyses on the original BIT-58 core between the surface 
and 126-m depth were carried out at Princeton University with a proce- 
dure modified from Bender et al.” and described by Higgins et al.. The 


standard deviation of *°Ar,,,, measured across 68 external standards 
is +0.014%o, corresponding to an age uncertainty of +210 kyr (10). For 
Allan Hills samples drilled in 2015-2016, *°Ar,,,, was measured at Scripps 
Institution of Oceanography with slight modification. First, sample size 
was increased from 500 to 800 g for better precision. Second, pumping 
time for evacuating the ice-bearing vessel before gas extraction was 
reduced to 15 min to prevent gas loss. Third, extracted and purified gases 
were measured ona newer mass spectrometer, which offers better signal 
sensitivity and stability. Finally, following published procedures*®, we 
measured the Xe/Kr ratios, expressed as 5Xe/Kr, an indicator of mean 
ocean temperature”, in the same aliquot of the extracted and getterred 
gases. A typical *°Ar,,, and 8Xe/Kr measurement takes ~3.5 h on the 
mass spectrometer. 

Two measurement campaigns at Scripps Institution of Oceanography 
were carried out and each was associated with slightly different ana- 
lytical procedures. In the first campaign, a fixed amount of reference 
gas was introduced into the bellow of the mass spectrometer, leaving 
different amount of gas in sample and reference sides. As the sample 
and reference gases deplete during analysis, the pressure and ion cur- 
rent fall more rapidly for the side containing less gas. As a result, Ar pm 
was observed to depend on the volume of the analyte relative to the 
reference. All measured *°Ar,,,, values were subsequently corrected for 
volume differences, using the voltage readings of the mass-40 (*°Ar) 
beam when the gases were introduced to the fully expanded bellow on 
the mass spectrometer. The magnitude of this correction is generally less 
than 0.02%, or 300 kyr, and the uncertainty in the correction is much 
smaller. Inthe second campaign, the amount of gas in the standard-side 
bellows was matched to that of the reference-side bellows. No volume 
correction was applied in the second campaign. Based on replicate 
analyses, the external reproducibility was +0.006%o (10) and +0.007%o 
(lo) for dry La Jolla air samples and Taylor Glacier blue-ice samples, 
respectively. The samples measured at Scripps therefore have an age 
uncertainty of +110 kyr (10) owing to the analytical uncertainty, or 10% 
of the sample age owing to the uncertainty in the *°Ar outgassing rate, 
whichever is greater. Individual sample uncertainties are reported along 
with the measured values. 

Xenon-to-krypton ratios (6Xe/Kr) were measured immediately after 
the *°Ar,mn Measurements on the same mass spectrometer by ‘peak 
hopping’ between “Kr and Xe. The 5Xe/Kr measurement takes about 
halfan hour to complete. The final reproducibility is +0.21%o and +0.27%o 
for dry LaJolla air samples and Taylor Glacier ice replicates, respectively. 
In order to calculate mean ocean temperature from 8Xe/Kr, which is nor- 
malized to LaJolla air, we need to make an assumption about the volume 
of the ocean”. For this purpose, we use a sensitivity of 1.1+ 0.2 °C/%o. We 
obtain this value from the observed last glacial maximum (LGM)-present 
8Xe/Kr change” with minor modifications. We adopted all the same sea- 
level-related corrections as in ref.” except one: the assumption that Kr 
and Xe saturation state was higher during the LGM than at present. This 
assumption lowers the sensitivity of mean ocean temperature (MOT) to 
6Xe/Kr from 1.2to 1.0 °C/%o (for example, a MOT change of 2.57 + 0.24 °C 
for a2.5%o change of 8Xe/Kr”). To account for this uncertainty we adopt 
an intermediate value of 1.1 and a larger error of +0.2. 


5D,,.and 5"°O,,,. Stable water isotope measurements (5'%0,.., 5D;-.) were 
taken from a slab of ice, 2.5-cm wide, cut (parallel to the vertical axis) 
from the outside of each ice core. Because some portions of the core 
were fractured, continuous sections were not available at all depths. 
These slabs were sub-sampled at 10-15-cm resolution at the Climate 
Change Institute, University of Maine. Each sub-sample was melted in 
asealed plastic bag at room temperature, vigorously shaken, and then 
decanted into an 18-ml plastic scintillation vial. All vials were refrozen 
and stored below -10 °C until the time of analysis. Water isotope samples 
were measured via cavity ring-down spectroscopy (CRDS) using a Picarro 
Model L2130-i Ultra High-Precision Isotopic Water Analyzer coupled 
with a High Precision Vaporizer and liquid auto sampler module. 6%0O,,, 
and &D,,. were measured simultaneously with an internal precision (20) 
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of 0.05% for 6'80,,., and +0.10%o for 6D,,.. The instrument drift adjust- 
ment and daily calibration were accomplished by periodically measuring 
three secondary standards: BBB (an average Maine freshwater fromthe 
Bear Brook watershed), ASS (Antarctic surface snow) and LAP (light 
Antarctic precipitation). These standards were calibrated against the 
internationally accepted water isotope standards: SMOW (standard 
mean ocean water), SLAP (standard light Antarctic precipitation) and 
GISP (Greenland ice sheet precipitation). 


5®N, 5°O.4m: and 5O,/N,. Analyses for the original BIT-58 O,/N,/Ar el- 
emental ratio, 80 of O,,and &°N of N, were carried out as described*””. 
Samples from ALHIC1502 and ALHIC1503 were measured using a slightly 
modified procedure. In brief, ~20 g of ice was placed into a glass flask 
chilled in a dry-ice-isopropanol cold bath (-78.5 °C) and the vessel was 
turbo-pumped for 3 min. The final pressure of the residual gas inside the 
vial was no higher than 1.5 mTorr. The pressure of the air trapped inice, 
when expanded to the same space, was ~3 Torr. The reduced pumping 
time thus introduces, at most, about 0.05% of modern atmosphere into 
the sample gas. 

After melting, the gas—meltwater mixture was equilibrated for at 
least 4h. The majority of the meltwater was drained and leftover water 
inside the glass vial refrozen at -30 °C. Once all meltwater had refrozen, 
headspace gases were collected by condensation at liquid helium tem- 
peratures. The sample was admitted to a mass spectrometer (Thermo 
Finnegan Delta Plus XP) for elemental and isotope ratio analysis after 
homogenizing at room temperature for 30 min. The reproducibility of 
8N, 60,/N, and 680, Values from ALHIC1502 and ALHIC1503 samples 
analysed using this procedure, calculated as pooled standard deviation, 
is +0.008%o, +3.164%o, and +0.024%o, respectively. 

Overall, the reduced pumping time leads to a reduction in gas loss 
and an improvement in precision. We note that the gas loss processes 
that alter 60,/N, did not appear to affect the isotopic composition of 
O, and N, in our samples (Extended Data Fig. 3). As a result, no gas loss 
correction was applied to the gas isotopes, in contrast to some previous 
studies in heavily fractured ice (for example, Siple Dome*’). All 6°N, 60,/ 
N, and 6'80,.m Values were normalized to modern atmosphere. 


CH,, CO, and 85°C of CO,. CH, concentrations were analysed at Oregon 
State University using a melt refreeze technique“. Samples (-60-70 g of 
ice) were trimmed, melted under vacuum, and refrozen at about -70 °C. 
CH, concentrations in released air were measured by gas chromatog- 
raphy and referenced to air standards calibrated by National Oceanic 
and Atmospheric Administration Global Monitoring Division (NOAA 
GMD) onthe NOAA04 scale. Precision was generally better than +4 ppb. 
Individual sample uncertainties, where available, are reported along 
with the measured values. 

CO, concentrations, referenced to air standards calibrated by NOAA 
GMD on the World Meteorological Organization (WMO) scale, were 
measured using a dry extraction (crushing) method”. Wherever pos- 
sible, ice samples were processed and analysed in replicate for each 
depth and results averaged to obtain final CO, concentrations. However, 
four pairs of replicates did not come from exactly the same depth and 
their reproducibility showed substantial deterioration beyond usual 
analytical uncertainties. In this case, we treated those four pairs of rep- 
licates as eight individual, unreplicated data points. For greenhouse gas 
concentrations, no gravitational fractionation correction was applied; 
the correction would be less than 0.8 ppb for CH, and 0.5 ppm for CO,. 

The measurement protocol for °C of CO, was as described” with 
improvements inthe efficiency of air extraction. Approximately 200 g ofice 
was dry-extracted with a custom-made ‘ice grater’. The reproducibility of 
the method, based on replicate analysis of Taylor Glacier blue-ice samples, 
is 0.02%o for 6°C-CO,. The method also provides a measurement of CO, 
concentrations. The uncertainties are +2 ppm (lo of the pooled standard 
deviation). The final °C values, normalized to the Vienna-Pee Dee Belem- 
nite (VPDB) standard, were corrected for blanks, N,O* mass interference 
from N,O, and gravitational fractionation using the nearest SN value. 


Age assignment of CO,/CH, and O,/N,/Ar samples 

The depths of samples analysed for CO,/CH, and O,/N,/Ar were differ- 
ent from depths of the *°Ar,,,, Samples that provided the chronology. 
CO,/CH, samples, O./N,/Ar samples, and ice samples were assigned ages by 
bracketing *°Ar,.n Values. Here we discuss two different models for assign- 
ing ages and showthat different criteria for age assignment lead to similar 
conclusions. For the purpose of investigating the relationship between 
greenhouse gases (CO, and CH,) and Antarctic climate (6D,,..), we binned 
samples into three units: MPT, 800-1,200 ka; pre-MPT, >1,200 ka; and young 
ice, <800 ka. We focus onthe robustness of ages for CO, and CH, samples. All 
CO, and CH, data discussed below exclude samples affected by respiration. 


Conservative approach. In the most conservative approach, an age 
was assigned only to samples bracketed by two *°Ar m ages lying with- 
ina single age bin. Otherwise, the sample was left ‘unassigned’. In this 
binning scenario, 19 CO, samples were classified as MPT, 15 as pre-MPT, 
and 34 as unassigned. We classified 28 and 19 CH, samples as MPT and 
unassigned, respectively. Eight CH, samples were binned into pre-MPT. 


Proximity approach. In this approach, each CO,/CH, datum was 
assigned the *°Ar,,,, age of the closest Ar isotope sample. In this age 
model, 23 and 37 CO, samples were binned into the MPT and pre-MPT 
units, respectively. For CH,, 31 samples were classified as MPT, and 16 
samples as pre-MPT. Finally, eight CO, and eight CH, samples were 
binned into the younger (<800 ka) age unit. 

Extended Data Fig. 4 summarizes the CO,-6D,,, and CH,-6D,.. 
relationships using the two age models. Our conclusions regarding the 
greenhouse gas-8D,,, relationship are not dependent on which model 
is used. As the proximity approach assigns an age to every depth, we 
adopted it for our interpretations. 


Assessing the fraction of climate variability recovered by 
discrete sampling 

When interpreting data from stratigraphically disturbed ice, the ques- 
tion arises as to what fraction of the amplitude of climate variability is 
accessed in the suite of samples analysed. Here we consider two factors: 
(1) the fraction of climate variability preserved in the ice cores compared 
to the true natural variability; and (2) to what extent discrete samples 
capture the variability that is preserved in the ice. Below we discuss 
these two factors separately. 


Fraction of the pre-MPT 5D,,, range preserved in the ice. We start 
by considering what fraction of the true glacial-interglacial variability 
in the 40k world is preserved in the ice. If we assume the absence of 
post-depositional alteration processes suchas respiration and diffusive 
smoothing, this fraction should be an intrinsic property of the ice itself 
and insensitive to the property being measured. Given the large number 
(>500) and the high spatial resolution (10 cm) of 6D,,. measurements, 
we used ôD; to evaluate the fraction of true 6D,,. variability preserved 
in the ice. Next, we needed to make an assumption of the true 6D,,, 
amplitude inthe 40k world. It is assumed that the true amplitude of 6D, 
variations in Allan Hills ice scales linearly with the amplitude of deep 
ocean temperature. This assumption is justified for two reasons. First, 
most of the bottom water formation today occurs in the southern ocean. 
Second, Antarctic temperature has been found to correlate tightly with 
mean ocean temperature during the LGM-Holocene transition” and 
with South Pacific bottom water temperature through the past 800 kyr“. 

For deep ocean temperature, we considered two scenarios. First, we 
assumed a constant partitioning between ocean temperature and the 
influence of seawater 80 on the benthic foram 680 record through- 
out the Pleistocene. This means that the amplitude of 6D,.. should 
simply scale with the benthic foram 680 values. The amplitudes of gla- 
cial cycles in the early and mid-Pleistocene (between marine isotope 
stages (MIS) 39 and 104), expressed in benthic foram 680, averages 
0.70%o, representing 37% of the change from MIS 6 to MIS 5e (1.90%0). 


The range of D,,. values between MIS 6 and MIS 5e observed in the Allan 
Hills ice from the 100k world was 80%o. Multiplying this 1OOk-world range 
by the scaling factor (37%) based on benthic foram 50 predicts the 
range of 6D,,, in 40k-world ice to be 30%o. By comparison, the range of 
40k-world 6D,,. observed in the Allan Hills ice is 32%o. Using this approach, 
we estimated that our pre-MPT ice samples represent -100% of the aver- 
age amplitude of glacial-interglacial 5D,.. cycles in the 40k world. 

In an alternative scenario, we considered deep ocean temperature 
inferred from a seawater 5'80 record constrained independently from 
records of carbonate microfossil 880 from the Mediterranean Sea*. 
Aslightly higher variability in the 40k world deep ocean temperature, 
which corresponds to 50% of the average variability in the 100k world, 
lowers our estimate of the recovered climate variability to -80%. Conse- 
quently, although other factors might influence the recovery of the range 
of climate variability, given these two approaches our best estimate is 
that the Allan Hills ice samples presented here preserve 90% of the true 
range of the climate variables in the 40k world. 


Fraction of the pre-MPT CO, and CH, range recovered by discrete 
samples. CO, and CH, concentrations were measured on a relatively 
small number of samples. It is therefore possible that we have not 
retrieved the full range of variability of these properties in the 40k world, 
evenif their true variability is well preserved in the ice cores. We evaluated 
this possibility by first generating a synthetic CO, series based upon the 
assumption that the variation of CO, in the 40k world scales linearly with 
benthic foram &'°0 values. A synthetic CO, time series between 1.2 and 2.0 
Ma was constructed on the basis of the LRO4 global benthic foram 680 
stack2. A linear regression between benthic 880 and atmospheric CO, 
between O and 800 ka resulted in correlation parameters that were used to 
create an artificial atmospheric CO, record further back in time (Extended 
Data Fig. 5a, b). After resampling to interpolate the synthetic data at a tem- 
poral resolution of 2.5 kyr, a Monte Carlo simulation randomly selected 
37 samples from the synthetic CO, record. The simulation was repeated 
10,000 times. We then calculated the observed range of the 37 random 
CO, samples and compared it to the ‘true’ range, yielding a ratio for each 
simulation run. Finally, the distribution of this ratio is shown as a series of 
histograms (Extended Data Fig. 5c-h). Note that the term of merit for each 
simulation is the ratio of the sampled range to the hypothetical true range. 

To account for the fact that interglacial accumulation rates are gener- 
ally higher than glacial rates*®, we multiplied the modelled occurrence of 
interglacial ice (CO, >250 ppm) by2, 4, 8, 16, and 32. When the presence 
of the interglacial ice is two times the glacial ice, 84 + 15% (95% confidence 
interval) of the total CO, range can be captured by 37 samples (Extended 
Data Fig. 5d). This ratio was derived from Vostok, an inland coring site 
of East Antarctica”. 

We note that 35 out of the 37 40k-world CO, samples involve replicates 
measured at the same depth. The assumption here is that replicates at 
the same depth have exactly the same age, which remains unevaluated 
in Allan Hills ice cores. It is possible that the ice flow and dip of ice layers 
make samples from the same depth different in age. If this is the case, 
each pair of replicates should be treated as two individual measure- 
ments, and the total number of samples would be 72. In this scenario, the 
recovered range is from 214 to 281 ppm, and 72 discrete samples capture 
88 + 13% of the CO, variability preserved in the ice. Treating replicates 
as individual samples does not significantly change the range (from 65 
to 67 ppm), but it increases the likelihood that we are recovering close 
to the full range of CO, in the ice (from 84% to 88%). 

Considering the discussion above, our best estimate is that the 37 
CO, samples cover 76 + 14% of the true CO, range in the 40k world 
(90% of true variability preserved in the ice multiplied by 84 + 15% that is 
captured by discrete sampling). Furthermore, we note that the 5D;,. 
values of the ice samples analysed for CO, as well as CH, span more than 
90% of the observed 40k-world 6D,,, range (Fig. 4). This number is higher 
than the estimate yielded by Monte Carlo simulations, meaning that the 
discrete samples are likely to capture more variability than probability 
alone would dictate. As a result, there is a greater chance that we have 


captured a significant fraction (>90%) of the glacial—-interglacial CO, 
range in the 40k world. 

For 6'°0,, Values the estimated recoverability should be no less than 
70%. First, 29 samples were analysed for 6'80,,,,. Second, the co-depth 
8D,.. values of 8%O,+m Samples span 97% of the 40k-world &D,,. (Extended 
Data Fig. 6a). 


Potential sample mixing on various length scales 

Mixing by molecular diffusion. The interval between 123 and 141m 
in ALHIC1503 shows unexpectedly low variability in the elemental and 
isotopic ratios of the major gases (Extended Data Fig. 7). Molecular 
diffusion in the firn column has been shown to reduce annual variability of 
water isotopes“ and gases”. In polycrystalline ice, diffusive mixing has also 
been invoked to explain the lack of sub-millennial 8D; variability preserved 
inthe EPICA Dome Cice core, dated -780 kyr old”. When annual layers are 
thin and temperatures are high, some gas records (for example, 60,/N,) 
can also undergo extensive diffusive smoothing and are expected to lose 
even the glacial-interglacial variability”. 

However, diffusive smoothing of climate records is unlikely to be 
significant at the Allan Hills BIA. First, current temperatures within the 
boreholes at the Allan Hills BIA are around -30 °C”. Rates of diffusion 
at the Allan Hills BIA should be substantially slower than at the warm 
basal temperatures (>-10 °C) expected at the bottom of thick (>3 km) 
polar ice sheets. Second, as is shown below, water isotopes (8D,,,) are 
expected to be much more sensitive to diffusive smoothing than gases 
(O, and N2). As a result, the presence of substantial sample-to-sample 
variability in 8D; values measured at 10-cm resolution indicates that the 
effects of diffusion on this and longer length scales are likely to be minor. 

Below we quantitatively estimate the characteristic diffusion time 
scale of water isotopes and of oxygen and nitrogen gases in the ice. 
Gas and water permeation coefficients in the ice have been estimated 
from molecular dynamics simulations°° and limited observations in 
some ‘natural experiments”. Values obtained from different methods 
can differ by more than one order of magnitude. In our calculation, the 
fastest permeation parameters are used whenever possible to provide 
the most rigorous test of diffusive smoothing. 

For the self-diffusion of water molecules in the ice, the diffusion time 
scale (Tice) depends solely on the diffusivity of water molecules in ice 
(Dice) and the length of interest (L; 0.1m): 
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Tice = D. (1) 
ice 


The diffusivity of ice is determined by™: 


Q. 
Dice = Dire * exp- re) (2) 


where D2, (9.1x10* m?s”) isa constant of diffusion for water molecules 
in ice, Qice (5.86 x 10*) mol") the activation energy of diffusion, R (8.314 
J mol” K”) the gas constant, and 7 the temperature. The characteris- 
tic time scale of self-diffusion in ice as a function of temperature is 
shown in Extended Data Fig. 8. Depending on temperature, it takes 10° 
to 10’ years for diffusive smoothing to smear the original signal on the 
length scale of 0.1m. 

Next, we consider the molecular diffusion of O, and N, in the ice. 
Below is a simple model of gas permeation witha large reservoir of gases 
(that is, bubbles) and a medium through which diffusion takes place 
(that is, ice): (1) First, gases in the reservoir dissolve into the ice. The 
amount of dissolved gas depends on ice volume, pressure, and solubility. 
(2) Next the gas will diffuse in the ice driven by the concentration gradient. 
This is assumed to be the rate-limiting step. (3) Diffusive gas exchange 
in the ice along the concentration gradient will immediately be com- 
pensated by the gas exchange between the ice phase and the bubble 
phase. (4) Eventually the concentration gradient in the bubbles will be 
balanced owing to diffusive flux. 
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Were there no gas reservoirs, the characteristic diffusion time scale 
of gas species m would simply follow the diffusivity of gas m in ice (Dm) 
and the length of interest (L): 


Tn= = (3) 


where L is 0.5 m because the spatial resolution of gas measurements 
is about 50 cm. 

Similar to the diffusivity of water isotopes, the diffusivity of gas m 
follows™: 


Dy = D°, x exp| “a 


Q. 
zi (4) 


where D? is a constant and Q,, cis the activation energy of diffusion. 
However, the presence of the reservoir would compensate for the 
diffusive flux. Here, we introduce the concept of partitioning function 
Z, the ratio of the number of gas molecules in the bubbles to the number 
of gases dissolved in the ice. 
E 


n= Zim * p (5) 
m 


By definition, Zof a given gas species m is 


n 
Zasa (6) 
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where Nn gas AN Np ice are the number of gas molecules in the gas phase 
and in the ice phase, respectively. 

Itis assumed that Nn gas > Nm ice- The number of gas m molecules in the 
gas phase can be inferred from the molar fraction of gas m (f,,) and the 
total number of gases in the ice (n,a), which has indeed been measured 
as the gas content of the ice (V), expressed in SI units m° kg: 


Mm gas = “gas * fy (7) 


Given the ideal gas law, equation (7) can be re-written as: 


PVmice 
RT? 


xf, (8) 


Nm, gas = 


where m,,, is the mass of the ice, p° is 101,325 Pa and 7° is 273 K. 

The number of molecules of gas m in the ice phase depends on the 
solubility of gas m (S„), the partial pressure of gas m (P„), and the number 
of water molecules (Nice): 


Nm ice = Sin x Pn x Nice (9) 
The solubility of gas m is governed by: 


Sm= 52x exo|- 9 (10) 


where s2 [Pa] is the dissolution constant, and Q,, the activation 
energy of dissolution. $5, and Sy, are 3.7 x 10 Pa“ and 4.510" Pa‘, 
respectively; Q), and Q, are 9,200 J mol” and 7,900 J mol", 
respectively®, ° í 

The partial pressure of gas m is the product of the bubble 
pressure, which is assumed to be the hydrostatic pressure of the 
depth in which bubbles are located, and the molar fraction of the gas 
minair: 


Pra= Pice %8 X 1) Xf, (11) 


where picis the density of ice (920 kg m°), g the gravitational accelera- 
tion constant (9.8 m s°), and A the depth. In our case, h is assumed to 
be a constant 150 m. 

The number of water molecules ni can be directly calculated from 
the mass of the ice: 


ice (12) 


where Mater iS the molecular weight of H,O (0.018 kg mol’). 
Therefore, the final expression for Z,,, is: 


Zn aro “Jn (13) 
mM s.x xoxhxf. x ise 
m Pice £ Ín Mwater 
10) 
= p Z Vx Myater 
m 0 (14) 
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Using equations (4) and (14) and relevant parameters, the diffusion time 
scale of O, and N, are computed as a function of absolute temperature 
(Extended Data Fig. 8). Unless at very cold temperatures (<230 K), water 
isotopes measured every 10 cm should be more susceptible to diffusive 
mixing than gases measured every 50 cm. As a result, the preservation 
of 8D,-. variability argues that flat 0,/N, and 6°0,.m profiles are not the 
result of diffusive smoothing. However, if the temperature is indeed that 
low, the characteristic time scale of gases becomes too large (>10’ yr) 
for substantial smoothing to occur in million-year-old ice. We did not 
calculate the permeation coefficients for CH, or CO, as these molecules 
have larger molecular diameters than O, and N,, and should diffuse at 
an even slower rate and havea longer t (ref. ”). 


Mixing within individual samples. In polar ice sheets the thickness of 
annualice layers declines with depth owing to compression and lateral 
flow. Thinning will reduce geochemical variability if a property changes 
over timescales shorter than that encompassed by the physical length of 
asingle sample. This problem is further amplified in stratigraphically dis- 
turbedice, as the orientation of the layering is uncertain. CO, and O,/N,/ 
Ar samples were cut from smaller lengths ofice than 8D,,. samples (which 
were 10-cm long). Therefore, assuming the layering is horizontal, CO, 
and O,/N,/Ar samples should average shorter time intervals than is the 
case for 5D,,.samples. On the other hand, CH, samples were cut in10-cm 
lengths. Because 6D,,, samples show large variability, we consider that 
gas properties other than *°Ar,,,, and 5Xe/Kr represent relatively short 
climate intervals without much mixing. The 6D,,, samples themselves 
may sometimes be mixed, and this possibility remains to be examined 
by the analysis of samples smaller than 10 cm. 


Flow-induced mixing. Physical mixing of ice of different ages can result 
from glacial flow. Mechanical mixing of this type may create an artifi- 
cial correlation between properties that are otherwise uncorrelated. 
Granted, many climate properties will co-vary in the absence of mixing. 
Therefore, property-property correlations do not necessarily indicate 
that mixing has occurred. 

Extended Data Fig. 9a shows the Pearson correlation coefficient 
(expressed as R°) between 8D, and d within intervals of varying lengths. 
In normally ordered ice cores, dis found to have a complex relation- 
ship with 6D,,.°°. We calculate R? as a function of the number of 5D,,. 
samples included. The section between 132 and 142 m shows a very strong 
correlation between 6D,,. and d, which raises the suspicion of 
mixing. 

We next examined the gas records between 132 and 140 m. The maxi- 
mum and minimum CO, values within this 10-m section are 253 and 
217 ppm, respectively. Because the gas content in the deep ALHIC1502 


and ALHIC1503 samples is very similar (-0.098 cm? g”), a simple linear 
mixing curve in the gas phase is expected. If we regard these two points as 
end members and plot 6D,,,. against CO,, the data points do not fall onto 
a mixing line (Extended Data Fig. 9b). Similarly, it is difficult to explain 
the 5D;,.-CH, and CO,-CH, plots by simple two-end-member mixing. 
In other words, a simple mixing scenario with two end members cannot 
explain our gas data, measured at every ~50 cm. We therefore consider 
the property-property plots (Fig. 4) to be unaffected by mechanical 
mixing onthe length scale of 50 cm, although mixing at a much smaller 
length scale is still possible. 


Effect of respiration on gas properties 

Respiration of detrital organic matter has been shown to lead to the 
in situ production of CO, in the bottom of polar ice sheets”, result- 
ing in elevated concentrations of greenhouse gases compared to the 
contemporaneous atmosphere. We observed substantially elevated 
CO, and CH, concentrations in the basal sections of ALHIC1503 and 
ALHIC1502, which we attribute to respiration and methanogenesis, 
respectively. Two lines of evidence support this hypothesis. First, there 
are extremely depleted (<-500%o) 60,,/N, ratios and high (>4%o) 68O tm 
values in samples near the bottom of ALHIC1503. For comparison, 
typical 5O,/N, values in Antarctic ice cores vary between -5 to -20%o* 
and6&80,,, ranges from -0.5 to +1.5%o” on glacial-interglacial timescales. 
Second, measured 8”C of CO, values in some deep samples from the 
ALHIC1502 and ALHIC1503 cores (Extended Data Fig. 10) are outside 
the range (-6 to -7%o) observed for the last 160 kyr”. This is consistent 
with contributions from respiratory CO, witha 86?C value of about —25%o. 
The deepest measured 6°C sample comes from 3 m above the bedrock 
in ALHIC1503 and has a 6°C value of -22.4%o. Assuming an initial 6°C 
of CO, value of -6.5%o and an 8"C of CO, value of -25%o from respired 
carbon, isotope mass balance suggests that 86% of the CO, inthis sample 
is derived from respiration. 

In four out of the nine 8°C samples, 6°C of CO, was compatible with 
the absence of respiration, assuming a cut-off 6°C value of -7%o. In the 
other five samples, which originated within 7 m of the bottom of the two 
cores, however, 86°C of CO, was lighter than -7%o and indicates respira- 
tion. We marked the shallowest 68°C sample with less than —7%o isotopic 
as a cut-off depth and rejected all CO, samples below it. Admittedly it is 
still possible that samples immediately above those cut-off depths might 
also be affected, but the highest pre-MPT CO, value is found around 
131m in ALHIC1503, a depth that is bracketed by two °C measurements 
falling between -6% and -7% and is considered respiration-free. As a 
result, even if certain samples above the threshold depths are contami- 
nated by respiration, our conclusions would not change. On the basis of 
these results, along with anomalous values of CH,, 50,/N,, and 880m 
inthe deeper sections, we excluded biogenic gas data below 185.950 m 
in ALHIC1502 and 139.625 m in ALHIC1503 from our evaluation. We still 
included 5D,,, from those sections in Fig. 2. 
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Extended Data Fig. 1| Satellite imagery of the Allan Hills study area. Right, 
WorldView03 colour pan-sharpened imagery (copyright 2011, DigitalGlobe, 
Inc.) of the mainice field, Allan Hills Blue Ice Area (in true colour mode). Inset, 
Antarctica with hill shading (source: Australian Antarctic Data Centre, Map 
13469; licensed under a Creative Commons Attribution 3.0 Unported License) 
on which our study area is marked by the black square. The black outcrop inthe 
main image is the Allan Hills nunatak. The red arrow marks the local iceflow. 
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Left, a magnified image of the drilling site from same source file. The imagery is 
enhanced (gamma-adjusted on ArcGIS 10) to highlight the colour contrast 
within the blue ice, which now has a greenish hue owing to the colour rendition. 
The brown lineations in the ice are exposed dust bands, providing a first-order 
tracer of surface ice stratigraphy. The locations of the cores reported in this 
work (see text) are marked with red circles. The location of the GPR profile in 
Extended Data Fig. 2 is shown asa yellow dashed line. 
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Extended Data Fig. 2 | GPR profile proceeding from the SW to the NE, 
crossing (within less than 5 m) ALHIC1502 and ALHIC1503. The exact transect 
location is shown in Extended Data Fig. 1. The location and depths of boreholes 
ALHIC1502 and ALHIC1503 (red bars) and ice drilled previously as BIT58 (grey 
bar) are indicated. The profile was collected with an 80-MHz MLF antennaina 
step-and-collect survey style with a step size of 25 cm and a stacking rate of 

64 scans. Standard post-processing steps were applied, including time-zero 
position correction, distance normalization, 50/110-MHz finite impulse 
response filter, background removal, and gain adjustment. Depth estimates 
from two-way travel time (TWTT) and migration use a radar travel velocity of 


0.165 m ns”. a, Non-migrated radargram showing dipping bed-parallel englacial 
stratigraphy and a strong dipping apparent bedrock reflection. Dashed blue line 
indicates the modelled location of the actual bedrock location as shown in b. 

b, Migrated radargram showing the ‘true’ location of the bedrock. Migration 
quality is poor at this location owing to the steeply dipping slope of the bedrock 
rise. Comparison between aand bsuggests that the correct bedrock reflector 
locationis translated downward and more steeply dipping than indicated inthe 
non-migrated data. ALHIC1502 and ALHIC1503 borehole depths independently 
verify this interpretation. 
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Extended Data Fig. 3 | Minimal impact of gas loss on oxygen and nitrogen 
isotopes. Pair differences of 850m (ASSO em; red) and SN (ASYN; blue) are 
plotted against A5O,/N, in ALHIC1502 and ALHIC1503 samples. The difference 
between two replicate samples is attributable to gas loss. The very weak 
dependence of A&80,,,, and ASN on ASO,/N, suggests that gas isotopes are 
relatively insensitive to gas loss fractionation in Allan Hills ice. 
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Extended Data Fig. 4 | Age assignment of CO, and CH, samples. Using the two 
age models discussed in the text— the conservative approach (top and bottom 
left) and the proximity approach (top and bottom right)—CO, and CH, are 
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plotted against 5D,... Both of the age models show reduced range inthe MPT and 
pre-MPT ice. We note, however, that the highest CO, and CH, values fall into the 
unassigned category in the strictest age model. 
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Extended Data Fig. 5 | Evaluating the fraction of CO, range captured by 37 time series (see Extended Data Fig. 5b) by 37 randomly selected samples, given 
discrete samples. a, Ice core CO, record between 0 and 800 ka versus LRO4 different preferential preservation of interglacial ice. Here the term ‘interglacial’ 
benthic foram "O stack synchronized onto the AICC2012 timescale. The is operationally defined as any sample with greater than 250 ppm CO.. Different 
result ofa linear regression is shown and used to calculate the synthetic CO, multipliers indicate the multiple occurrence of the ‘interglacial’ ice to simulate 
record between 1.2 and 2.0 Ma. b, A synthetic CO, record between 1.2and2.0Ma varying degrees of ice preservation biases. This analysis shows that when the 
reconstructed from LRO4 680 and the regression parameters calculated ina. presence of interglacial ice is no more than eight times greater than that of 


The range of CO, in this synthetic time series is 213-269 ppm.c-h,Fractionofthe glacial ice (c-f), we could expect 37 random samples to be likely to capture 
recovered CO, variability (observed range/true range) of the synthetic CO, 66-98% (95% confidence interval) of the CO, range in the ice. 
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Extended Data Fig. 6 | Evaluating the fraction of 5O m and 5Xe/Kr range 
captured by discrete samples. a, b, Allan Hills 5°O,,,, (a) and 6Xe/Kr (b) are 
potted against the 5D,.. of the same depth, colour coded according to their age 
units. Vertical dashed lines represent the range of 5D,.. observed in the 
40k-world ice (-284 to -316%o). Notably, the range of 6D,,. associated with nine 
5Xe/Kr values from the 40k world is only about 40% of the entire 6D,.. rangein 
our 40k ice samples (-288 to -300%o); low 8D;ce values (less than -300%o0) are 
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missing from the 6Xe/Kr samples. Thus, we donot expect the nine Xe/Kr samples 
to fully capture the range of Xe/Kr ratios in the 40k-world ice. On the other hand, 
the co-depth §D,,, values of 8°O,. Samples occupy 97% of the total range, 
implying that 29 640, Samples are likely to cover 70% of the variability 
preserved in the 40k-world ice, barring any diffusive smoothing of the “O,,,, 
records. 
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Extended Data Fig. 7 | Gas and ice properties in the interval between 123 and all measurements: 0.010%o for 6°N, 3.063% for 60,/N, and 0.026%o for 5O tm- 
143 min ALHIC1503. a, 8°N; b, 50,/N,;¢, 5O m; d, 5D;,.. The error bars Note that the range of &°0,,,,, is less than 0.2%o. By comparison, glacial-to- 
associated with gas properties represent the pooled standard deviation (10) of interglacial 6O m variability in the 100-kyr climate cycles is about 1.5%0”". 
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Extended Data Fig. 8 | Evaluating the effect of diffusive mixing onice and gas 
properties. a, Partitioning function Z (ratio of the number gas molecules in the 
gas phase to that in the ice phase) of O, (blue) and N, (red) asa function of 
temperature. b, Characteristic time scale of the self-diffusion of water molecules 
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inice (black), of O, permeation in ice (blue), and of N, permeation in ice (red), 
plotted as a function of temperature. Note that for water molecules the length 
scale of diffusion (L) is 0.1m, whereas for N, and O, L =0.5 m, reflective of their 
respective sampling resolutions. 
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Extended Data Fig. 9 | Evaluating the possibility of flow-induced mixing in 
Allan Hills ice cores. a, A contour map shows the correlation coefficient (R°) 
between 6D,,,and din core ALHIC1503 as a function of the number of adjacent 
stable water isotope data points included in the calculations, from5 to 85 witha 
step of 10. The sampling resolution is approximately 10 cm. The high correlation 
coefficient (>0.7) between 132 and 142 m raises the suspicion that mixing has 
taken place in this interval. However, the possibility of mixing is not supported 
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by CO, and CH,. b-d, Cross-plots of 6D;.e-CO, (b), 5D;-e—CH, (c), and CO,-CH, (d) 
in intervals between 132 and 140 m in ALHIC1503. The black solid lines are mixing 
lines based on two end members at 135.30 m and at 139.76 m, where the 
minimum and maximum CO, values are measured, respectively. Most of the 
points do not fall onto the mixing line, implying that two-end-member mixing 
alone cannot explain the high correlation between 5D,,.and d. 
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plottedinred, and along with 8°C, plotted in blue) are also shown for 
comparison. Note the marked scale difference for CO, and 6"C between 


ALHIC1503 and ALHIC1502. 


Extended Data Fig. 10 | Respiration in basal ice revealed by 6°C of CO, in 
ALHIC1503 and ALHIC1502. a, ALHIC1503; b, ALHIC1502. 5D,., (top), CO, 
(middle), and 86°C of CO, (bottom) in Allan Hills ice cores are plotted against 
depth. 5D,..and CO, (both measured independently in smaller samples, 


